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Abstract 

A  finite  difference  method  for  solving  the  time- 
dependent  Navler-Stokes  equations  for  an  Incompressible 
fluid  Is  Introduced.   This  method  uses  the  primitive 
variables,  i.e.  the  velocities  and  the  pressure,  and  is 
equally  applicable  to  problems  in  two  and  three  space 
dimensions.   Test  problems  are  solved,  and  an  application 
to  a  three  dimensional  convection  problem  is  presented. 


Introduction,   The  equations  of  motion  of  an  Incompressible 
fluid  are 

^,u.  +  u.S.u.  = h.p   +   vV^u.  +  E.     (V^  =  y        S^) 

S  .u  .  =  0 

J  J 

where  u.  are  the  velocity  components,  p  Is  the  pressure,  p 
Is  the  density,  E.  are  the  components  of  the  external  forces 
per  unit  mass,  v  Is  the  coefficient  of  kinematic  viscosity, 
t  Is  the  time,  and  the  Indices  l,j  refer  to  the  space 
coordinates  x . ,  x  .,  1,  j  =  1,2,3«   S-  denotes  differentiation 
with  respect  to  x. ,  and  h.    differentiation  with  respect  to 
the  time  t.   The  summation  convention  is  used  in  writing 


the  equations. 

We  write 


^    d'=^  ^  d 


where  U  is  a  reference  velocity,  and  d  a  reference  length. 
We  then  drop  the  primes.   The  equations  become 

d.u.  4-  Ru.S.u.  =  -  S.p  +  V^u.  +  E.         (1) 

S.u  -  0  (2) 

J   J 

where  R  =  Ud/v  is  the  Reynolds  number.  It  is  our  purpose 


to  present  a  procedure  for  solving  these  equations  by  finite 
differences  In  some  bounded  domain  X^  In  either  two  or  three 
space  variables.   The  distinguishing  feature  of  this  method 
Is  the  use  of  equations  (1)  and  (2)  rather  than  higher  order 
derived  equations.   This  makes  It  possible  to  satisfy  all 
boundary  conditions  and  to  achieve  adequate  computational 
efficiency  even  in  problems  involving  three  space  dimensions 
and  time.   The  author  is  not  aware  of  any  other  method  for 
which  these  claims  ca,n  be  made. 


Principle  of  the  method.   Let  u. ,  p  denote  the  solution  of 
(1)  -  (2)  as  well  as  its  discrete  approximation,  and  let 
Du  =  0  be  a  difference  approximation  to  ^ -u  .  =  0.   In 
general  we  shall  allow  Du  to  take  different  forms  in  the 
interior  of  the  domain  &  and  on  its  boundary;  at  the  boundary 
we  may  wish  to  use  higher  order  one-sided  differences  so 

as  to  preserve  accuracy  of  some  given  order.   It  is  assumed 

n   n 
that  at  time  t  =  nAt  velocity  and  pressure  fields  u^,  p 

are  given  satisfying  Du"^  =  0.   The  task  at  hand  is  to 

evaluate  u?"^"^,  p^"^''"  using  equation  (1),  so  that  Du"^   =  0. 

Auxiliary  fields  u?"^,  are  first  computed  through 

uj^  =  u^;  +  At  F.u  (3) 


where  F.u  approximates 


-  Ru.a.u.  +  V^u.  +  E.     '  V^  =  3    a^ 
J  J  1       11         "  ^ —  J 


F.u  may  depend  on  u. ,  u.   ,  and  intermediate  fields, 
say  ui^,  u^*.   In  the  evaluation  of  u.    equation  (2)  and 
the  pressure  term  In  equation  (1)  are  not  taken  into  account. 
An  iteration  procedure  is  now  Introduced  for  the  purpose  of 
finding  u.    inside  B-   and  p    in  >^  and  on  its  boundary, 
by  setting 

n+1,1    n  /i,  X 

P      =  P  (4a) 

n+l,m+l  aujx         .^^m  ,  /),,_\ 

Uj^      '  =   u^        -  AtG^p  m   >  1  (4b) 

^n+l,m+l        _n+l,m        ,  ^^  n+l,m+l  ^  -,  , ,,    . 

p         '  =  p         '       -   ADu         '  m   >  1  (4c) 

where  A  is  a  parameter,  the  quantities  u?   '"^  and  p^"*"-*-'"^ 

are  successive  approximations  to  u^   ,  p"^   ,  and  g'?p  Is  a 

^    ,  .     f,     n+l,m    -,   n+l,m+l   ,  .  , 

function  of  p    '      and  p    '     which  converges  to  a 

difference  form  of  S.p"^^^  as  jp'^+l'^i+l  _  p'^+^^'^l  tends  to 

zero.   The  form  of  <a%   is  crucial  for  the  convergence  of 

the  iterations  and  for  the  accuracy  of  the  over-all 

scheme;  it  will  be  specified  below. 

When  for  some  S,   and  a  small  predetermined  constant  e 

max  |p^+l'^+l  -  p^+l'^l  <  e  (5) 

we  set 

^n+1  ^  ^n+1,^+1 
1      1 

pn+1  ^  pn+1,^+1 


If  the  velocity  component  u.  Is  prescribed  at  the  boundary, 
u^  '^       In  (4b)  and  (4c)  Is  replaced  by  the  given  value  of 
u.   .   The  Iterations  (4)  Insure  that  equation  (1),  Including 
the  pressure  term.  Is  satisfied  Inside  X^,    and  equation  (2) 
Is  satisfied  In  ^  and  on  Its  boundary. 

The  proposed  method  Is  essentially  a  generalization  of 
the  artificial  compressibility  method.  Introduced  In  [  1  ] 
for  the  purpose  of  finding  steady  solutions  of  equations 
(1)  -  (2).   This  can  be  seen  In  the  following  way:   Let  us 
drop  the  requirement  (5)  and  set  £  =  1   for  all  n.   t  then 
does  not  represent  real  time,  since  equation  (2)  Is  not 
satisfied  and  we  are  not  following  the  evolution  In  real 
time  of  the  solutions  of  (1)  -  (2).   We  choose  for  G.p 
an  expression  which  approximates  B.p,  and  set  5  =  At/A. 
Remembering  the  definition  of  u.   ,  we  see  that  equations  (4) 
represent  one  step  In  an  artificial  time  for  a  system  approxi- 
mating the  equations 


2 
3,u.  +  Ru  .S  .u.  =  -  ^.p  +  V  u.  +  E. 
ti      Jjl      i       li 


S,  p  +  ^  S  .u  . 


(6) 


This  is  the  auxiliary  system  used  in  [  1  ]  to  find  steady 
solutions  of  equations  (1)  -  (2).   If,  as  the  artificial 
time  is  advanced,  the  solution  of  (6)  tends  to  a  steady 
limit,  i.e.  one  which  does  not  depend  on  t,  then  this  solution 
is  a  steady  solution  of  equations  (1)  and  (2). 


Comparison  with  the  method  of  Harlow  and  Welch.   Before 
Introducing  specific  scheme  for  evaluating  u.    and  a 
specific  representation  for  Du  and  G^^p,  we  shall  compare  our 
method  to  the  only  other  method  known  to  the  author  for 
solving  the  incompressible  Navier- Stokes  equations  using 
the  primitive  variables.   That  method  is  due  to  Harlow  and 
Welch  [  2  ],  who  applied  it  to  two  dimensional  problems.   The 
purpose  of  the  comparison  is  to  explain  what  the  iterations 
(4)  accomplish.   It  should  be  noted  that  in  an  Important 
respect  the  method  of  Harlow  and  Welch  is  more  ambitious 
then  ours,  having  been  applied  to  flows  with  a 
free  surface,  while  we  have  thus  far  attempted  to  solve 
only  problems  for  a  confined  fluid.   Disregarding  for  the 
moment  the  details  of  the  difference  scheme  used  by 
Harlow  and  Welch,  and  using  notations  similar  to  those 
introduced  in  the  preceding  section,  we  can  summarize  their 
method  as  follows: 

Let  Du  =  0  approximate  ^  .u  .  =  0,  and  G.p  approximate 
S.p.   It  is  assumed  that  at  time  t  -   nAt  velocity  fields 
u.  are  given,  satisfying  Du  =  0.   Equation  (2)  can  be 
approximated  by 

u"?^^  =  uj  +  AtLuJ  -  AtQ^u'^  -  AtQ^p'^  +  AtE^       (7) 

where  Lu  approximates  V  u,  and  Q.u  approximates  S.u.u.. 

J-  J  -^  J 

Performing  the  operation  D  on  (7),  and  assuming 
Du"+1  ^  0 


we  obtain  the  following  equation 

l'p^  =  -  §r  +  ^Lu""  -  DQu^  +  DE.  (8) 

I  2 

where  L  p  =  DGp  approximates  V  p.   This  equation  is 

of  course  a  difference  analogue  of  the  equation 

V^p  =  -  a.S,u.u  +  S  E  (9) 

which  can  be  obtained  from  equation  (1)  by  taking  its 
divergence.   Equation  (8)  is  solved  by  iteration;  the  boundary 
conditions  are  obtained  by  applying  equation  (1)  to  the 

fluid  at  the  boundary.   The  pressure  field  p   thus  obtained 

n+1 
Is  used  in  equation  (7)  to  evaluate  u.   ;  we  then  have 

T.  n+1   ^ 
Du    =0. 

Our  method  naturally  bears  a  resemblance  to  the  method 

of  Harlow  and  Welch.   By  substituting  equation  (4b)  into 

(4c),  which  is  possible  except  near  a  boundary,  we  obtain 

the  equation 

pn+l,m+l  _  pn+l,m  ^  _  ^^^aux  ^  ^^  ^^mp 

which,  in  view  of  the  definitions  of  D,  G^  and  u^   ,  is  an 
iteration  procedure  for  solving  an  analogue  of  equation  (8). 
In  this  sense  our  method  is  related  to  Harlow  and  Welch's 
like  a  predictor-corrector  method  to  a  predictor  method; 

whereas  Harlow  and  Welch  first  determine  p   so  that  Du    =0 

.  J  ,     n     ^   n+1   n+1     J  j-L, 
we  make  a  guess  at  the  values  at  u   ,  p    ,  and  then 

n+1 
correct  them  until  the  constraint  Du    =  0  is  satisfied. 
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The  advantage  of  the  predlctor-corrsctor  approach  is  that 
it  allows  us  to  evaluate  u  .   by  Implicit  schemes,  and 
thus  take  much  larger  time  steps;  this  is  of  course 
particularly  useful  if  three  dimensional  problems  are  to 
be  solved. 

Harlow  and  Welch  retain  the  terms  -rr-  Du  ,  DLu  =  LDu 
in  equation  (8),  although  it  assumed  that  Du  =  0.   This 
allows  them  to  reduce  the  stringency  of  the  convergence 
criterion  in  the  iterative  solution  of  (8),  and  thus  to 
reduce  the  computational  effort.   If  we  consider  the 
definition  of  u.   ,  we  see  that  this  feature  of  Harlow  and 
Welch's  method  is  preserved  in  ours  in  a  natural  way; 
we  shall  see  in  a  later  section  that  through  an  appropriate 
choice  of  G.p  we  can  Introduce  further  time-saving  devices. 

It  is  believed  that  the  main  advantage  of  our  formulation 
over  that  of  Harlow  and  Welch  lies  in  the  following:  the  system 
consisting  of  equation  (1)  and  (2)  and  the  one  consisting 
of  (1)  and  (9)  are  equivalent  only  if  equation  (2)  is 
satisfied  at  all  boundaries.   In  order  to  ensure  that  (2) 
is  satisfied  Harlow  and  Welch  introduce  artificial  reflection 
principles  at  the  boundaries;  in  effect,  this  means  that  the 
assigned  boundary  data  are  not  imposed  on  equation  (2);  the 
solution  thus  obtained  has  therefore  an  uncertain  relation 
to  that  of  the  problem  to  be  solved.   In  our  formulation,  no 
such  problem  arises;  equations  (1)  and  (2)  are  consistently 
approximated  everywhere,  and  the  correct  boundary  data  used. 
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Harlow  and  Welch  implement  their  method  with  a  difference 
scheme  which  contains  many  attractive  features,  and  which  -. 
could  be  used  with  our  formulation.   The  schemes  we  are 
about  to  introduce  appear  to  be  much  more  efficient  than 
theirs,  but  suitable  mainly  for  a  problem  where  the  domain 
/d'  has  a  simple  shape,  and  where  the  boundary  data  are  smooth. 


o  1  i"Y 

Evaluation  of  u. .   We  shall  now  present  schemes  for 

evaluating  u.   ,  defined  by  (3)- 

In  [  1  ]  we  used  a  combined  Duf ort-Prankel  and  leap-frog 

scheme,  in  which  time  and  the  first  space  derivatives  were 

approximated  by  centered  differences,  and  a  second  derivative 

such  as 

S-,u 

was  replaced  by 

This  scheme  is  appropriate  only  when  an  asymptotic  steady 
solution  Is  sought.   It  is  however  inacceptably  inaccurate 
when  real  time  dependence  is  studied,  unless  At  is  inadmlssibly 
small.   Our  reason  for  discussing  this  scheme  here  is  the 
following:   The  Duf ort-Prankel  scheme  is  explicit  and 
unconditionally  stable;  it  is  therefore  a  natural  scheme  to 
use  when  the  non-linear  terms  in  (1)  are  differenced  in 


"conservation-law"  form.  I.e.  In  the  form 


S  .  (  u.  u  . ) 
J   1  J 


rather  than 


u  .^  .u.  . 
J  J  1 

We  found  that  in  problems  in  which  the  viscosity  is  not 
small,  it  is  preferable  to  use  "non-conservative"  difference 
schemes  for  the  non-linear  terms,  and  avoid  the  Dufort- 
Frankel  scheme. 

We  shall  instead  be  using  one  of  the  following  variants 
of  the  alternating-direction  implicit  method: 

A.   In  two-dimensional  problems  we  shall  use  a 
Peaceman-Rachford  scheme,  as  proposed  by  Wilkes  [ 3  ]  in 
a  different  context.   This  takes  the  form 

*      _  n        p  At   n     ,  *  *       . 

^l(q,r)  -  ''i(q,r)  "  ^  5a7^  ^1  (q,r )  ^  ^i{q+l,  r )  "  ^i{q-l,r)^ 

n  At   n     /  n  n       « 

"  ^  WKI^   "2(q,r)^''i(q,r+l)  "  ^i(q,r-l)^ 

At   /  *  *  o   *  \ 

^   i^  ^^i{q+l,r)  +  ^i{q-l,r)  "  ^^i(q,r)) 

At   /  n        ,   n  ^,  n     » 

i^  ^^l(q,r+l)  +  "i(q,r-l)  "  2^i(q,r)^ 

+  —  E^ 


(10a) 
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^l(q,r)  =  ^^Kq^r)  "  ^  WK^   "l(q,r )  ^^l{q+l,r )  "  "i(q-l,r)) 

T3  At   *     /  aux        aiix     x 
"  ^  T^   ^2(q,r)^"i(q,r+l)  "  ^i(q,r-l)^ 

At   /  *  *  *      v 

2^^^^2  ^  i(q+l,r)     i(q-l,r)     i(q,r)' 

At   /  aiix         aux        ^  aux    . 

2A?  ^""K^^^+D   ""iCq.^-D  "   i(q.r)^ 

+  —  E^ 

where  u.  are  aiixlliary  fields,  and  u.  ,    s    =   u.  (qAxj_,rAxo)  • 
As  usual,  the  one -dimensional  systems  of  algebraic  equations 
can  be  solved  by  Gaussian  elimination. 

B.  In  two  dimensional  and  three  dimensional  problems 
we  shall  use  a  variant  of  the  alternating  direction  method 
analyzed  by  Samarskii  in  [  4  ].   This  takes  the  form 

*        _  n  P  At   n       /  *         _   * 

^i(q,r,s)  -  ^i(q,r,s)  "  "  2A^  ^l(q,r,  s  )  ^^i  (q+1,  r ,  s  )   ''l(q-l,r,s) 

At   /  ■*  * 

+  -^   ("i(q+l,r,s)  +  ^i(q-l,r,s) 


-  ^^l(q,r,s)) 


^l(q,r,s)  =   ''i{q,r,s)  "  ^  ^^   ''2(q,r,  s)  ^ ''i(q,r+l,  s  )"  ''i{q,r-l,  s  )  ^ 

At    **  ** 

+  ^  ("i{q,r+l,s)  +  ^i{q,r-l,s) 


-  2u. 


i{q,r,s) 
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aux  ,  **  R     At       **  '  r.i^^  1,3^  ^ 

'^i(q,r,s)  i(q,r,s)  2Ax        ^iQjr-jS)'    i(q,r,s+l)         i(q,r,s-l)' 

,     At     /    aux  ,      aixx  ^   auj<  ^ 

^^2    ^    i(q,r,s+l)  i(q,r,s-l)  i(q,r,E)' 

+  AtE. ,  , 

i(q,r,s) 


^l(q,r,s)    ■=  ^i(qAx^,rAx2,sAx^) 
^l(q,r,s)    ^   E.(qAx^,rAx2,sAx3) 


u.  ,  u.   are  airxlllary  fields.   These  equations  can  be  written 
In  the  symbolic  form 

(I  -  AtQ^)u*  =  u"^ 

(I  -  AtQ2)u''''  =  u*  (11) 

(I  -  AtQ^ju    =  u 

where  I  Is  the  identity  operator,  and  Q.  involves  differentia- 
tions with  respect  to  the  variable  x.   only.   It  can  be 
verified  that  this  scheme  is  accurate  to  0(At).   If  R  /  0, 
schemes  (10)  and  (11)  are  accurate  to  the  same  order. 
Scheme  (11)  is  stable  in  three  dimensional  problems,  while 
we  do  not  know  of  a  simple  extension  of  (10)  to  the  three 
dimensional  case.   Scheme  (11)  has  two  useful  properties: 
It  requires  fewer  arithmetic  operations  per  time  step  than 
scheme  (10),  and,  because  of  the  simple  structure  of  the 
right  hand  sides,  the  intermediate  fields  u. ,  u.   do  not 
have  to  be  stored  separately. 
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If  the  scheme  (11)  is  to  be  used  in  a  problem  In  which 
the  velocities  u^  are  prescribed  at  the  boundary,  values 
of  u.,  u.  ,  u    at  the  boundary  have  to  be  provided  so  that' 


the  operators  (I  -  AtQ. )  can  be  inverted.   The  scheme  is 
accurate  to  0(At),  and  a  variety  of  choices  of  u. ,  u.  , 


u?^^  can  be  made  without  detracting  from  that  accuracy;  for 
example,  we  can  set 

aux  *  4(-*  n 

u.         =  u.    =   u.      =  u. 

1  111 


or 


aux  *  **  n+1 

u.         =   u.    =  u.      =  u. 
1  111 


However,  in  the  iteration  (4)  we  need  the  values  of  the 
partial  differences  of  u?^^,  and  it  can  be  seen  that  with  a 
careless  choice  of  u.,  u.  ,  u.    at  the  wall  the  error  in  Du 

would  be  0(At/Ax).   To  ensure  that  Du  ^^  is  accurate  to  0{At) 

*   **   aux 
we  use  boundary  values  of  u. ,  u.  ,  u.    consistent  with  the 

scheme  (11),  for  example 

u*  =  uf  1  -  AtQ^u^;-^^  -  AtQ^uJ+l  ^  ^^Q^pn 


(12) 


**    n+1    ^,  „   n+1  ,  .^.n   ^^ 
u.   =  u.    -  AtQ,u.    +  AtG.p 


,  aux    n+1  ,  .,  „  n 
u^   =  Uj_   +  AtG^p 
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where  G.p  approximates  S.p  to  0(At).   We  assume  throughout 
that 

At  =  0(Ax?)  . 


By  Inspecting  equations  (11)  we  can  see  that  no 
reasonable  results  can  be  expected  unless 

RAt  «  1  ;  (13) 

this  Is  an  eminently  reasonable  condition;  In  view  of  the 
definitions  of  R  and  of  the  dlmenslonless  time,  It  merely 
states  that  during  a  given  time  step,  a  typical  fluid  blob 
covers  a  distance  which  Is  small  compared  to  a  typical 
dimension  of  the  problem.   The  Inequality  (13)  suggests 
that  as  R  Is  Increased,  the  time  step  At  has  to  be  decreased. 

It  should  be  noted  that  our  non-dlmenslonallzatlon  Is 
not  the  standard  one;  our  At  corresponds  to  a  real  time 
interval  R  times  longer  than  the  one  which  corresponds  to  At 
In  the  standard  non-dlmenslonallzatlon 

f  =  (§)t  . 
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As  an  Illustration,  we  applied  these  three  schemes  to  a 
simple  parabolic  problem.   We  solved  the  heat  equation 


h^u   =  i^l   +  hl)u 


in  the  two  dimensional  square  ^  :  0  <_  x-,  <_  tt,  ^  ;^  ^^  ^  Tr> 
with  the  initial  data 

u(0)  -   sin  X-,  sin  x^ 

and  the  boundary  condition  u  =  0  on  the  boundary  of  ^. 
For  this  problem,  u    =  u   ;  all  three  schemes  considered: 
Dufort-Frankel,  (10)  and  (11)  are  stable  for  all  values  of 
At.   The  Dufort-Prankel  scheme  however  is  consistent  only 
when  At  =   o(Ax. ) . 

The  exact  solution  of  the  above  problem  is 

-2t 
u(t)  =  e    sin  x-,  sin  Xp 

In  Table  I  we  display  the  maxima  of  the  differences  between 
the  exact  solution  and  the  solutions  obtained  using  respectively 
the  Dufort-Frankel  scheme  and  schemes  (10)  and  (11).   The 
Dufort-Frankel  scheme  requires  two  levels  of  initial  data; 
we  set 

u(At)  =  e~    sin  x,  sin  Xp 

in  accordance  with  the  exact  solution. 

We  chose  Ax-,  =  AXp  =  Ax  =  Tr/l9;  At/Ax   =1.   n  is  the 
number  of  time  steps.   We  see  that  for  this  value  of  At 
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Table  I:   Solutions  of  the  heat  equation 


n 

Dufort-Frankel 

A 

B 

1 

0.00011 

0.00080 

2 

0.0^57 

0.00021 

0.00152 

5 

0.0408 

0.00030 

0.00216 

4 

0.0314 

0.00038 

0.00273 

5 

0.0208 

0.00045 

0.00323 

6 

0.0111 

0.00051 

0.00368 

7 

0.0028 

0.00057 

0.00406 

8 

0.0039 

0.00062 

0.00440 

9 

0.0093 

0.00066 

0.00469 

10 

0.0134 

0.00069 

0.00493 

20 

0.0204 

0.00080 

0.00574 
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the  Dufort-Frankel  scheme  Is  not  accurate,  and  displays 
oscillations  which  can  be  explained  by  the  presence  of  the 
term 

fe^  ^t^ 

in  the  expression  of  the  truncation  error  for  that  scheme. 
Scheme  (10)  is  more  accurate  than  scheme  (11);  we  shall 
see  that  this  advantage  is  not  as  notable  when  the  non-linear 
terms  are  present  in  the  equations. 


The  Dufort-Frankel  Scheme  and  successive  point  over-relaxation. 
In  order  to  describe  our  iteration  procedure  for  determining 
the  pressure  and  satisfying  the  equation  of  continuity,  we 
need  a  few  facts  concerning  the  Dufort-Frankel  scheme  for 
the  heat  equation  and  its  relation  to  the  relaxation  method 
for  solving  the  Laplace  equation. 
Consider  the  equation 

-  V^u  =  f,      V^   ^  hf   +  hf  (14) 

in  some  nice  domain  ^,    say  a  rectangle,   u  is  assumed  known 
on  the  boundary  of  /^.   We  approximate  this  equation  by 

-  Lu  =  f  (15) 

where  L  is  the  usual  five  point  approximation  to  the 
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Laplaclan,  and  u  and  f  are  now  m-component  vectors,   m  Is 
the  number  of  Internal  nodes  of  the  resulting  difference 
scheme.   For  the  sake  of  simplicity  we  assume  that  the 
mesh  spacings  in  the  x-,  and  Xp  directions  are  equal. 
Ax,  =  AXp  =  Ax;  this  implies  no  essential  restriction.   The 
operator  -  L  is  represented  by  an  m  x  m  matrix  A  which  is 
positive  definite,  symmetric,  and  has  positive  diagonal 
elements. 

We  write 

A  =  A'  -  E  -  E' 

where  E,  E'  are  respectively  strictly  upper  and  lower 
triangular  matrices,  and  A'  is  diagonal.   The  convergent 
relaxation  iteration  scheme  for  solving  (I5)  is  defined  by 

(A'  -  aaE)u'^"^-'-  =  {(l-^)A'  +  a)E')u^  +  cof         (I6) 

(see  e.g.  [5]).   cd  is  the  relaxation  factor,  0  <  03  <  2, 
and  the  u^  are  the  successive  Iterates.   The  evaluation 
of  the  optimal  relaxation  factor  <^q„^  depends  on  the  fact 
that  A  satisfies  "Young's  condition  (A)",  i.e.  that 
there  exists  a  permutation  matrix  P  such  that 


P~-'-AP  =  A  -  N  (17) 


where  A  is  diagonal,  and  N  has  the  normal  form 

Vg'  oJ 
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The  zero  submatrlces  being  square.   Under  this  condition, 
"^opt  ^^^  ^^   readily  determined  (see  e.g.  [  5  ],  Chapter  9). 

The  matrix  A  depends  on  the  order  in  which  the  components 
of  u    are  computed  from  u"^.   Changing  that  order  is 
equivalent  to  transforming  A  into  P~  AP,  where  P  is  a 
permutation  matrix. 

We  now  consider  the  solution  of  (14)  to  be  the 
asymptotic  steady  solution  of 

S^u  =  V^u  +  f  (18) 

and  approximate  the  latter  equation  by  the  Dufort-Frankel 
scheme 

n+1    n-1   2At  /  n     ,  n      n      n       n+1   n-1 « 


^q,r  =   u(qAx-j^,rAX2,nAT) 


or,  grouping  terms 


Ax"^   ^'^        Ax"^  ^'^ 


-  2  5  ("^+l,r-q-l,rK,rrtK,r-l)  ^  ^^Tf 


Since  u    does  not  appear  in  (I9),  the  calculation 
q^  ■'■ 

separates  into  two  independent  calculations  on  intertwined 
meshes,  one  of  which  can  be  omitted.   When  this  is  done, 
we  can  write 
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U^        =1     pn+l  /        ^^  ^^^   "^  components) 

Vu 


If  we   then  write 

Q    At 


"        2 

144 

At 

03   =        ^^  ,  ^  (20) 


we  see  that  the  Iteration  (19)  reduces  to  an  Iteration  of 
the  form  (l6)  where  the  new  components  of  U    are  calculated 
in  an  order  such  that  A  has  the  normal  form  (l?)-   The 
Dufort-Frankel  scheme  appears  therefore  to  be  a  particular 
ordering  of  the  overrelaxation  method  whose  existence  is 
equivalent  to  Young's  condition  (A). 

The  best  value  of  At,  At   , ,  can  be  determined  from 
CD   .  and  relation  (20).   We  find  that  At   .  =  0(Ax),  and 
therefore  for  At  =  At   ,  the  Dufort-Frankel  scheme  approxi- 
mates, not  equation  (l8),  but  rather  the  equation 


V  =  ^'"  -  2(|l''  S?"  • 


This  is  the  equation  which  Garabedian  in  [  6  ]  used  to 

estimate  co  ^,  .   It  can  be  used  here  to  estimate  At   .  . 
opt  opt 

These  remarks  obviously  generalize  to  problems 
where  Ax-,  ^   AXp  or  where  there  are  more  than  two  space 
variables. 
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^    J  ^    .  .     n+1   n+1 
The  Iteration  procedure  for  determining  u. ,  p 

For  the  sake  of  clarity  we  shall  assume  in  this  section  that 

the  domain  /^Is  two  dimensional  and  rectangular,  and  that 

the  velocities  are  prescribed  at  the  boundary.   Extension 

of  the  procedure  to  three  dimensional  problems  Is  Immediate, 

and  extension  to  problems  with  other  tjrpes  of  boundary 

conditions  often  possible.   Stress-free  boundaries  and 

periodicity  conditions  in  particular  offer  no  difficulty. 

Let  3   denote  the  boundary  of  ^   and  C  the  set  of  mesh 

nodes  with  a  neighbor  in  ^.   In  ^  -  i5  we  approximate  the 

equation  of  continuity  by  centered  differences,  i.e.  we  set 

^"^  ^   2Ai^  ("l(q+l,r)  "  ^l(q-l,r)^ 

+  2^  (U2(q,r+1)  "  ^2{q,r-l))  =  ° 

At  the  points  of  «  we  use  second  order  one-sided 
differences,  so  that  Du  is  accurate  to  0(Ax  )  everywhere. 
Consider  the  boundary  line  Xp  =  0,  represented  by  j  =  1 
(Fig.  1).   We  have  on  that  line 

^^  "  A%  f^2(q,2)  -  ^2(q,l)  "  \    ("2(q,5)  "  ^2{q,l)^^ 

,   1  (22) 


(21) 


with  similar  expressions  at  the  other  boundaries. 

n-f-1  ' 
Let  (q,r)  be  a  node  in  ^  -  ;^  -  C,    knowing  u.   ' 

1  =  1,2,  and  p   '  ,  we  shall  evaluate  simultaneously 
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n+l,m+l  n+l,m+l      ^    -,      n+l,m    ,  _^     /    .  ,> 

l(q±l,r)^    ''2(q,r±l)    ^^^   Pq,r         (Pq,r=  pC^Ax^,  rAx^ ) )  ,    using 


u 

the  formulae 


^n+l,m+l    n,m       ^„  n+l,m+l 

Pq,r     =  Pq;r  -  ^^^  (23a) 

(D  given  by  (21) ) 

n+l,m+l      _     aujc                      At      ,    n+l,m        1  ,    n+l,m+l          n+l,m,.  ,„,,  , 

^l(q+l,r)    -   ^l(q+l,r)    "   2A5q;   ^V2,r   "   2  ^Pq,r             "  Pq,r      ^^  ^^^b) 

n+l,m+l      _     aux                      At      ,1    ,    n+l,m+l    ,    ^n+l,m«      ^n+l,m>  , ^^    . 

^l(q-l,r)    -   ^l(q-l,r)    "   2A57   ^2    (Pq,r  +  Pq,r      )"   Pq-2;r)  ^^^c) 

n+l,m+l      _     aux                      At      ,    n+l,ni        1  ,    n+l,m+l   ^     n+l,m^,  ,„,,, 

2(q,r+l)    -   ^2(q,r+l)    "   2A3^   (Pq,r+2    "  2  (Pq,r             +  Pq,r      ^^  ^^^d) 


u 


n+l,m+l      _     aux  At      ,1    ,    n+l,m+l    ,      n+l,mv        n+l,m,       ,„,    , 

^2(q,r-l)    -   ^2(q,r-l)    "  2A3^   ^2    (Pq,r  +  Pq,r      )"  Pq,rl2^      ^^^e) 

(See  Fig.    2). 

These  equations  specify  O^p   for  use  in  equation  (4b). 

In  (^   and  Q  these  formulae  have  to  be  modified.   Consider 
again  the  boundary  line  x^  =  0  (Pig.  1).   Assume  the  velocities 
are  prescribed  at  the  boundary,  i.e.  n}t     ,«  are  given, 
1  =  1,  2.   There  are  several  ways  of  including  that  information 
in  the  iterations  (4).   The  consistent  way  (in  the  sense 
of  equations  (12))  would  be  to  set 


and 


aux       n+1     ■  A-i-n  ^ 
u.  /    .  =  u. /   T \  +  AtG. p 
i(q>l)     i(q,l)      1^ 


n+l,m+l  _   n+1 
^i(q,l)   -  ^i(q,l) 
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This  is  the  better  approach,  which  we  intend  to  use  in  future 
work.    However,  for  the  sake  of  simplicity,  we  chose  here 
an  inconsistent  way  of  treating  the  boimdary,  i.e.  we  set 


i(q,l)     i(q,l)     i(q,l)      i(q,l)  ^ 

It  can  easily  be  verified  that  this  does  not  affect  the 
values  of  ul^  ;  it  does  however  introduce  an  additional 
error  of  0(At)  into  the  computed  pressure  field. 

Equations  (23)  and  their  modifications  near  the  boundary 
can  be  solved  for  p^+l^ni+l^  yielding,  for  (q,r)  in  ^-   Jg' -  d: 

n+l.m+1    /-,  i-lr/-,  \  n+l,m   ,T^  aux 

Pq  r     =  (1  +  a^^  +  02)   [  (1  -  a-j^  -  a2)p   '   -  ADu 

(25a) 
,  n+l,m  ,   n+l,mx  ,  ^  ,  n+l,m  ,  ^n+l,mv-, 
+  °l(Pq+2;r  +  Pq-S^p)  +  °^2(Pq,r+2  +  Pq,r-2)^  ' 


for  (q, r)  in  C 

Pq,2     =  (1  +  a-^  +  2  ^2^    [(1  -  a-L  -  ^  0'2^Pq,2    "  ^^"^ 

,    n+l,m  ,   n+l,mv  .  „  ^n+l,m-| 
+  °l(Pq+2;2  ^  V2;2)  +°'2Pq,4   ^ 

and  for  (q,r)  on  "23  : 

pn+l,m+l  ^  (,  ^  a,)-l[(l  -  a,)p^;i'^  -  >.Du' 


aux 


(25b) 


aux 


,  _   ,  n+l,m   1  ,_n+l,m    n+l,mwn 
+  2o^2(Pq,3   "  ^  ^Pq,4    "  Pq,2   ^^^ 


(25c) 


etc.   In  (25a)  and  (25b)  Du  is  given  by  (21),  and  in 

(25c)  by  (22).   The  several  representations  of  Du  have  been 

used  to  derive  these  expressions,  and  we  have  written 


25 


a.  =  ^     1  =  1,  2  . 

These  Iterations  are  to  be  performed  until  for  some  2 

max  I  n+l,i+l    n+l.-^i  < 
q,r  '^q,r       ^q,r   I  - 

for  a  predetermined  e. 

Formula  (25a)  Is  seen  to  be  a  Duf ort-Frankel  relaxation 

scheme  for  the  solution  of  an  analogue  of  equation  (8); 

this  of  course  was  to  be  expected.   The  at  of  the  preceding 

section  is  replaced  here  by  XAt/2j  corresponding  to  at  ^ 

(or  (JO  ^  )  we  find  A   ,  .   If  P  were  known  on  ^  or  <S  and  the 
^     opt  opt      ^ 

Iterations  (25a)  stood  alone,  convergence  of  the  iteration 
for  all  A  >  0  would  follow  from  the  argument  of  the 
preceding  section,  and  A  =  ^^nt  ^^^^"^  lead  to  the  fastest 
convergence.   Although  no  proof  is  offered,  the  numerical 
evidence  leads  us  to  state  that  the  whole  iteration  system — 
equations  (25a),  (25b),  (25c) — converges  for  all  A  >■  0,  and 
converges  fastest  when  A  =  A   , .   None  of  the  boundary 
instabilities  reported  by  users  of  the  stream  function 
vorticity  approach  has  been  observed. 

In  (25a)  the  Laplacian  of  p  is  evaluated  by  a  difference 
formula  using  a  stencil  whose  nodes  are  separated  by  2Ax-,, 
2AXp.   This  results  from  the  use  of  centered  differences 
for  approximating  both  S  .u .  and  S.p.  As   a  consequence,  the 
pressure  iterations  split  into  two  calculations  on  intertwined 


24 


meshes,  coupled  at  the  boundary.   The  most  efficient  orderings 
for  performing  the  Iterations  (25)  a^^e  such  that  the 
resulting  over-all  scheme  Is  a  Duf ort-Frankel  scheme  for 
each  one  of  the  Intertwined  meshes.   This  Involves  no 
particular  difficulty;  a  possible  ordering  for  a  rectangular 

grid  is  shown  in  Fig.  3. 

n+1 
The  new  velocities  u.   ,  1  =  1,  2,  are  to  be  evaluated 

using  (25a),  (23b),  (23c),  and  (23d).   This  has  to  be  done 

only  after  the  p    '   have  converged;  there  Is  no  need  to 

evaluate  and  store  the  intermediate  fields  u.   '    .A 

1 

saving  in  computing  time  can  be  made  by  evaluating  Du 
once  and  for  all  at  the  beginning  of  the  iteration. 
An  advantage  of  our  formulation  appears  here:   The  Iterates 
p    '   in  (25a),  (25b),  (25c)  converge  to  the  solution  of 
a  certain  set  of  equations.   Suppose  we  were  to  solve  this 

set  of  equations  by  some  other  iteration  procedure, 

^    .     u   ^  .    .^   ^    n+l,i+l   n+1,^ 

stopping  when  two  successive  iterates  p       ,    V 

differ  by  less  than  a  small  number  e;  we  could  then  use 

the  latest  iterate,  p    '    to  evaluate  u.   ,  though 

formulae  such  as 

,n+l  _  au:!C      n+1,^+1 
^i   -  ^1   ~  ^iP 

p      n+l 
where  G.p  approximates  S.p.   If  At  =  0(Ax  ),  Du    would 

fail  to  vanish  by  an  amount  0(e).   On  the  other  hand,  using 

formulae  (23),  when  p^+^>^+^   and  p"^^-*-'-^  differ  by  less 

than  e,  Du    =  0(e/A);  it  can  be  seen  that  A   ,  =  0(Ax~  ), 
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therefore  Du"^   =  O(eAx);  a  gain  in  accuracy  appears,  which 
can  be  used  to  relax  the  convergence  criterion  for  the 
Iterations. 


Solution  of  a  test  problem.  The  proposed  method  was  first 
applied  to  a  simple  two  dimensional  test  problem,  used  as 
a  test  problem  by  Pearson  [  7  ]  for  a  vortlclty-stream 
function  method.  J9'  Is  the  square  0  <  x^  <_  ir,  1  =  1,  2; 
E-,  =  Ep  =  0;  the  boundary  data  are 

-2t 
u,  =  -  cos  X-,  sm  Xp  e 

-2t 
Up  =   sm  X-,  cos  Xp  e 

The  Initial  data  are 

U-,  =  -  cos  x^  sin  Xp 

Up  =   sin  x-j^  cos  Xp 

The  exact  solution  of  the  problem  Is 

-2t 
.-2t 

-4t 


u,  =  -  cos  X-,  sm  Xp  e 


Up  =   sm  X-,  cos  Xp  e 


p  =  -  R  ^  (cos  2x,  +  cos  2Xp)( 
where  R  Is  the  Reynolds  number.   We  assume  that 


Ax^  =  AXg  =  Ax 


A   ^  Is  then  found  to  be 
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2 
•X     =      2Ax^ 
opt   At  sln(2Ax) 

In  tables  II,  III,  IV,  V  and  VI  we  display  the  results 
of  some  sample  calculations,   n  is  the  number  of  time  steps; 
e(u.),  i  =  1,  2,  are  the  maxima  over  -^  of  the  difference 
between  the  exact  and  the  computed  solutions  u. .   It  is  not 
clear  how  the  error  in  the  pressure  is  to  be  represented; 
p   Is  defined  at  a  time  intermediate  between  (n-l)At  and 
nAt;  it  is  proportional  to  R  in  our  non-dimensional - 
ization.   There  are  errors  In  p  due  to  the  use  of  the  inconsistent 
boundary  data  (24),  and  due  also  to  the  fact,  discussed  at 
the  end  of  the  preceding  section,  that  the  iterations  can 
be  stopped  before  the  p^'"^  have  truly  converged.   e(p)  in 
the  tables  represents  the  maximum  over  the  grid  of  the  differ- 
ences between  the  exact  pressure  at  time  nAt  and  the  computed 
p^,  divided  by  R;  it  is  given  mainly  for  the  sake  of 
completeness.   The  accuracy  of  the  scheme  is  to  be  Judged 
by  the  smallness  of  e{u.).   i  is  the  number  of  iterations; 
it  is  to  be  noted  that  the  first  iteration  has  to  be 
performed  in  order  to  approximate  equation  (1).  "Scheme  A" 
means  formulae  (10)  were  used  to  evaluate  u   ,  "Scheme  B" 
means  formulae  (11)  were  used.   e  is  the  convergence 
criterion. 
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Table   II 
Scheme   A;    Ax    =  7r/39;    At    =   2Ax^    =   0.01597;    e   =  Ax^;    R   =   1 


n 

( 

5(1 

a,) 

e(i 

a^) 

e(p) 

i 

1 

2.8 

X 

10-^ 

2.6 

X 

10-^ 

0.0243 

1 

2 

2.7 

X 

10-^ 

2.0 

X 

10-^ 

0.0136 

7 

3 

1.5 

X 

10-^ 

1.3 

X 

10-^ 

0.0069 

4 

4 

1.8 

X 

10-^ 

1.9 

X 

10-^ 

0.0145 

4 

5 

1-3 

X 

10-^ 

1.7 

X 

10-^ 

0.0089 

5 

6    . 

lO 

X 

10-^ 

1.8 

X 

10-^ 

0.0116 

4 

7 

1.6 

X 

10-^ 

1.9 

X 

10-^ 

0.0144 

4 

9 

1.4 

X 

10-^ 

1.7 

X 

10-^ 

0.0147 

4 

10 

1.3 

X 

10-^ 

1.6 

X 

10-^ 

0.0156 

4 

20 

1.8 

X 

10-^ 

2.3 

X 

10-^ 

0.0241 

4 
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Table   III 
Scheme   A;   Ax   =  7r/39;   At   =  2Ax^   =  0.01397;    e   =  Ax^;    R  =  1 


n 

e(u^) 

e(u2) 

e(p) 

& 

1 

8.5 

X 

10-5 

3.8 

X 

10-5 

0.0059 

10 

2 

1.0 

X 

10-^ 

5.7 

X 

10-5 

0.0067 

10 

3 

1.0 

X 

10-^ 

7.0 

X 

10-5 

0.0068 

10 

4 

1.0 

X 

10-^ 

7.8 

X 

10-5 

0.0068 

10 

5 

1.0 

X 

10-^ 

8.3 

X 

10-5 

0.0069 

10 

6 

9.7 

X 

10-5 

8.6 

X 

10-5 

0.0070 

10 

7 

9.4 

X 

10-5 

8.7 

X 

10-5 

0.0071 

10 

8 

9.0 

X 

10-5 

8.7 

X 

10-5 

0.0073 

10 

9 

8.7 

X 

10-5 

8.7 

X 

10-5 

0.0077 

10 

10 

8.5 

X 

10-5 

8.5 

X 

10-5 

0.0082 

10 

20 

1.0 

X 

10-^ 

1.0 

X 

10-^ 

0.0216 

9 
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Table   IV 
Scheme   A;   Ax   =  Tr/j>9;    At   =  Ax^   =   0.00648;    e   =  Ax^;    R   =  20 


a 

e(u^) 

e{u^) 

e{p) 

& 

1 

2.5  X  10"^ 

2.7  X  10"^ 

0.0167 

15 

5 

3.7  X  10"^ 

3.5  X  10"^ 

0.0175 

10 

5 

5.5  X  10"^ 

4.5  X  10"^ 

0.0183 

10 

7 

6.8  X  10"^ 

5.2  X  10~^ 

0.0189 

10 

9 

7.5  X  10"^ 

5.4  X  10"^ 

0.0193 

10 

:0 

7.0  X  10"^ 

4.7  X  10"^ 

0.0204 

9 
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Tatle  V 
Scheme  B;   Ax   =  Tr/39;   At   =  Ax2   =  0.00648;    e   =  Ax2;    R  =  20 


n                        e{u^)                            e{u^)  e(p)  Z 

1  3.9  X  10--^  4.4  X  10-5  0.0404  l6 

3  5.9  X  10-5  g^Q  X  10-5  0.0466  11 

5  8.5  X  10"^  6.7  X  10-5  0.0505  10 


e(u^) 

3.9   X 

10-5 

5.9   X 

10-5 

8.5   X 

10-5 

1.0  X 

10-2 

1.1    X 

10-2 

1.0  X 

10-2 

7       1.0  X  10-^       7.4  X  10"5      0.0551       10 
9       1.1  X  10""       7.9  X  10"5      0.0599       10 


20       1.0  X  lO"'^       7.8  X  10-5      0.0839       10 
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Table  VI 
Scheme  B;   Ax   =  ir/39;   At    =  i  Ax^   =   0.00324;    e   =  Ax^;    R   =   20 


n 


1 

e(u-^) 

e(u2) 

e(p) 

i 

1 

1.1    X 

10-5 

1.2   X    10-5 

0.0217 

15 

3 

1.9  X 

10-5 

2.1   X    10-5 

0.0234 

9 

5 

2.5  X 

10-5 

2.8   X    10-5 

0.0242 

9 

7 

3.3  X 

10-5 

3.2   X    10-5 

0.0249 

9 

9 

4.0  X 

10-5 

3.5  X   10-5 

0.0253 

8 

:0 

5.8  X 

10-5 

3.9  X   10-5 

0.0258 

8 
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Tables  II  and  III  describe  computations  which  differ 

2 
only  in  the  value  of  e.   They  show  that  e  =  Ax  is  an 

adequate  convergence  criterion.   Tables  IV  and  V  indicate 

that  reliable  results  can  be  obtained  even  when  the  Reynolds 

number  is  large  for  the  grid  employed;  when  R  =  20,  Ax  =  ^  , 

we  have 

R  2:  1-5  Ax~ 

The  errors  are  of  the  order  of  1  '^  .   Tables  IV  and  V  show 
that  scheme  B  does  not  introduce  unduly  large  errors  into 
the  calculations;  table  VI  shows  that,  as  expected,  a 
decrease  in  At  decreases  the  error  (but  increases  the 
amount  of  computational  labor). 


Application  to  thermal  convection.   Suppose  a  plane  layer  of 

fluid,  of  thickness  d  and  infinite  lateral  extent,  is  heated 

from  below.   The  lower  boundary  x^  =  0  is  maintained  at 

a  temperature  T  ,  the  upper  boundary  x^  =  d  at  a  temperature 

Tn  <  T  .   The  warmer  fluid  at  the  bottom  expands  and  tends 
1    o 

to  move  upward;  this  motion  is  inhibited  by  the  viscous 
stresses. 

The  equations  describing  the  fluid  motion  are,  in  the 
Boussinesq  approximation  (see  e.g.  [8  ],  Chapter  2) 
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-1         p     ' 
h.u.    +   u.S.u.  =   -  -^   S.p  +  vV  u.  -  g(l  -   (T  -  T  ))5. 

^,  T  +  u  .a  .T  =  kV^T  ^  .u  .  =  0 

t      J  J  J  J 

where  T  Is  the  temperature,  k  the  coefficient  of  thermal 
conductivity,  a  the  coefficient  of  thermal  expansion,  and 
5.  the  components  of  the  unit  vector  pointing  upwards. 
We  write 

d         '     ^-^1  v2 

-1  =  (v)-i       T   -r^        ^'  =  (J)* 

o  1 


1   d  2     (T.-T  )dgx 


'    X. 


and  drop  the  primes.   The  equations  now  are 

S,u.  +  u.S.u.  =  -  S.p  +  V^u.  +  5-  (T-l)5. 
a^T  +  u.a  .T  =   ^  V^T 

S  .u  .  =0 

J  J 

where  R*  =  agd-^(T  -  T-,  )/kv  Is  the  Raylelgh  number,  and 

0  =   v/k  the  Prandtl  niomber.   The  rigid  boundaries  are  now 
situated  at  x^  =  0  and  x^  =  1,  where  It  Is  assumed  that  u.  =  0, 

1  =  1,  2,  3. 

-X-      ^ 

It  Is  known  that  for  R  <  R  ,  the  state  of  rest  Is 

c' 

*         /- 
stable  and  no  steady  convection  can  arise,  with  R  =  1707-7d2 

(see  [9])-   When  R  =  R  ,  steady  Infinitesimal  convection 


3^ 


can  first  appear,  and  the  field  quantities  are  given  by 


u^  =  CW(x^)(t) 

C 

u.  = 


^  =  ^  W(x^)B^<t)    1  =  1,  2 


T  =  CT(x^)<t) 

where  ^   =  <t)(x-,,Xp)  determines  the  horizontal  planform 
of  the  motion  and  satisfies 

{^l   +  hi  +  a^)^  =   0  , 

W(x^),T(x^)  are  fully  determined  functions  of  x^,  a  =  3-117j 
and  C  is  a  small  but  undetermined  amplitude,   (see  [8]). 

In  two  dimensional  motion  u-,  =  0  and  the  motion  does 
not  depend  on  x, .   ¥e  then  have 

(j)  =  cos  aXp 

The  motion  is  periodic  in  x^  with  period  2Tr/a . 

The  Nusselt  number  Nu  is  defined  as  the  ratio  of  the 
total  heat  transfer  to  the  heat  transfer  which  would  have 
occured  in  no  convection  were  present.   For  R  ^  ^^ f 
Nu  =  1.   In  our  dlmensionless  variables 

27r/a 


Nu  =  ^  /      (au^T  -  a^T)dX2 


In  [1]  and  [10]  we  studied  the  dependence  of  Nu  on  R*  and 
o  for  steady  convection,  assuming  that  the  motion  remains 
periodic  in  Xp  with  period  2Tr/a.   There  seems  to  be  no  point 
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in  reproducing  these  calculations  with  the  present  method; 
some  values  of  Nu  were  recomputed  and  are  in  good  agreement 
with  those  obtained  previously. 

It  appears  of  interest  to  determine  the  onset  of  instabil- 
ity by  finite  differences.   We  assume  that  the  flow  is 
periodic  in  x   with  period  2Tr/a,  and  we  define  r  =  R  /R  . 
We  find  that  for  r  =  1,  Nu  oscillates  between  1  and  1.03 
for  a  very  long  time;  for  r  =  O.99,  Nu  =  1.0000;  i.e.  the 
computed  and  the  exact  values  of  R  agree  to  within  less 

than  1  ^  .   These  results  were  obtained  with  AXp  =  2'rr/28a, 

/  2  2 

Ax^  =   1/27j  e  =  Ax„ ,  At  =  3Ax^.   The  corresponding  A   ,  is 
P  '-J  opt 

2     2 
2,  2AXp  sin  (ttAx^) 

A_.  = ,   p  =  1 ^    ^ 


'°P'  '  (At  ^  A|.)/i_p2  '     "       (A.|  +  Ax2) 
Ax^   AXg 

In  three  dimensional  convection  problems  not  only  the 
amplitude  of  the  motions  is  to  be  determined  but  also  their 
spatial  configuration.   For  R  =  R  ^   can  be  any  periodic 
solution  of 

{'b\  +  'b\^  ^^)k  =  0       a  =  3-117 

It  is  reasonable  to  assume  that  the  cell  patterns  are  made 
up  of  polygons  whose  union  covers  the  (x^jXp)  plane. 
Possible  cell  shapes  are  hexagons,  squares,  and  rolls 
(i.e.  two  dimensional  convection  cells).   It  is  known 
(see  e.g.  [11])  that  the  equations  of  motion  admit  steady 
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solutions  corresponding  to  more  than  one  kind  of  cell. 

However,  only  cellular  structures  which  are  stable  with 

respect  to  small  perturbations  can  persist  in  natiire  or 

be  exhibited  by  our  method.   It  has  been  shown,  numerically 

by  the  author  [10],  experimentally  by  Koschmieder  [12]  and 

Rossby  [13],  theoretically,  in  the  case  of  infinite  cr  and 

small  perturbations,  by  Busse  [14],  that  for  R  /R  <  10  the 

stable  cellular  mode   is  a  roll  with  a  wave  number  in  a 

certain  range.   We  shall  now  demonstrate  numerically  the 

impermanence  of  hexagonal  convection  and  the  emergence  of  a 

roll. 

We  choose  R/R  =r=2,  cr  =  i.   We  assume  the  motion 
'  c        ' 

to  be  periodic  in  the  x-,  and  Xp  directions,  with  periods 
respectively  47r//3a  and  47r/a  (the  first  period  is  apparently 
in  the  stable  range  of  periods  as  predicted  by  Busse  [14]). 
These  are  the  periods  of  the  hexagonal  cells, which  could 
arise  when  R  =  R  .   The  state  of  rest  is  perturbed  by  adding 
to  the  temperature  in  the  plane  x.,  =  Ax^  a  multiple  of  the 
function  <t)(x,  ,x_)  which  corresponds  to  a  hexagonal  cell,  and 
adding  a  small  constant  to  the  temperature  on  the  line 
X,  =  ^  (47r/v/Ja),  Xp  =  •^  (47r/a).   We  then  follow  the  evolution 
of  the  convection  in  time.   We  use  a  net  of  24  x  24  x  25,  i.e. 


Ax-^  =  (4Tr/Aa)/24  , 
AXp  -  (47r/a)/24 
Ax^  =  1/24 
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2 

we  also  set  e  =  Ax- 


At  =  3Ax^ 


We  then  find 

4 


A 


°P*    /  At  ,  At  ,  At  ,  r~2~ 
( — ^  +  — 2    — 2^"^   P 
Ax^   AXp   Ax^ 

2   2    2 
Ax-,Ax„  sin  ttAx^ 

where      ^  -  i    o  -^   ^ 2 

r   "~  P    P         ^    P         P    P 

(Ax^Ax^  +  Ax^Ax^  +  Ax|Ax^) 

To  visualize  the  convection  pattern,  we  proceed  as 
follows:   We  consider  the  plane  x^  =  17Ax  ;  if  u^,    -,o\  >  0 
we  print  an  *,  if  u^^     -,  p  \  ^  0  we  print  a  0. 

The  evolution  of  the  convection  is  shown  in  Fig.  4a,  4b, 
4c,  4d,  4e  and  4f.   The  hexagonal  pattern  introduced  into 
the  cell  is  not  preserved.   The  system  first  attempts  to 
achieve  equilibrium  as  a  roll  with  period  47r/a  (Pig.  4c), 
then  as  a  roll  with  period  STr/ZJa  (Fig.  4d),  and  finally 
becomes  a  roll  with  period  47rX/5"a  (Fig.  4f). 

With  different  initial  perturbations  the  evolution  of 
the  system  is  different,  the  final  result  is  the  same.   The 
calculation  was  not  pursued  until  steady  convection  had  been 
established,  because  that  would  have  been  excessively  time 
consuming.   It  is  known  from  previous  studies  that  steady 
rolls  can  be  obtained.   (See  e.g.  [1]),  and  that  the  mesh 
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used  here  provides  an  adequate  representation. 


Conclusion  and  applications.   The  Benard  convection  problem 
is  not  considered  to  be  an  easy  problem  to  solve  numerically 
even  In  the  two  dimensional  case  (see  e.g.  [I5]).   The  fact 
that  with  our  method  reliable  time  dependent  results  can 
be  obtained  even  In  three  space  dimensions  Indicates  that  the 
Navler-Stokes  equations  do  Indeed  lend  themselves  to  numerical 
solution.   A  niimber  of  applications  to  convection  problems, 
with  or  without  rotation,  can  be  contemplated;  In  particular. 
It  appears  to  be  of  Interest  to  study  systematically  the 
stability  of  Benard  convection  cells  when  cr  ^  00  ,  and  when 
the  perturbations  have  a  finite  amplitude.   The  existence 
of  finite  amplitude  Instability  of  the  cells  Is  suggested 
by  the  numerical  results  reported  in  [10] • 

Other  applications  should  include  the  study  of  the 
finite  amplitude  instability  of  Polseuille  flow,  the  stability 
of  Couette  flow,  and  similar  problems. 

Another  type  of  application,  presently  being  carried 
out  by  the  author,  includes  the  finite  difference  solution 
of  the  Green-Taylor  problem  [16]  and  related  problems  involv- 
ing simplified  representations  of  turbulence. 
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Figure  I.       Mesh  Near  a  Boundary 
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Figure2.        Iteration  Scheme 


hi 


The  domain  is  sweeped  in  the  order  AB,CD,EF,  GH,  IJ,K 


Figure  3.     An  Ordering  for  the   Iteration  Scheme 
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Pig.  k.      Evolution  of  a  Convection  Cell 
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4a.   After  1  step 
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4b.   After  10  steus 
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4c.   After  125  steps 
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4e.   After  325  steps 
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